/**
 * @file NewtonInterp.H
 * @brief 一个广义Newton（即包括Hermite）插值的实现类，
 * 可以根据输入的插值条件输出差分表，插值多项式;
 * 可以在已有的基础上新增插值条件，合并构建新的差分表与插值多项式;
 * 可以根据输入的插值条件与目标点直接给出插值近似值。
 * @author htli (3180102114@zju.edu.cn)
 * @version 1.0
 * @date 2021-10-24
 * 
 * @copyright Copyright (c) 2021  linhuo2020
 * 
 */
#ifndef _PROJECT2_NEWTONINTERPHT_H_
#define _PROJECT2_NEWTONINTERPHT_H_

#include <iostream>
#include <vector>

#ifndef _TYPEDEF_REAL_
#define _TYPEDEF_REAL_
///根据需要更改 double 为 long double 或 float
typedef double Real;
typedef std::vector<Real> RealVect;
typedef std::vector<RealVect> RealBox;
#else
//do nothing
#endif

#include "InterpConditions.H"
#include "Polynomial.H"
#include "TriangleTable.H"

const int fac(const int x)
{
    int f;
    if (x == 0 || x == 1)
        f = 1;
    else
        f = fac(x - 1) * x;
    return f;
}

class NewtonInterp
{
protected:
    TriangleTable<Real> tableOfDividedDiffs_;
    Polynomial<Real> interPoly_;
    InterpConditions IntpConds_;

public:
    NewtonInterp(){};
    NewtonInterp(const InterpConditions &_IntpConds) { this->ResetIntpConds(_IntpConds); };
    ~NewtonInterp(){};
    void ResetIntpConds(const InterpConditions &_IntpConds) { IntpConds_ = _IntpConds; };
    const InterpConditions getIntpConds() const { return IntpConds_; };
    ///在不求插值多项式的前提下求得某点处的估值（必须是狭义Newton型插值，否则报错并返回-1）
    static Real NevilleAitken(const Real &_x, const InterpConditions &_IntpConds);
    ///根据插值条件创建差分表，插值多项式
    void QuickInterp();
    void GenerateNewtonPoly();
    ///在已有的基础上增加新的插值点及对应条件，合并构建新的差分表与插值多项式
    void AddIntpSiteThenInterp(const RealVect &_SiteInfo);
    void MergeNewtonPoly(const int &_newAddedRowNum);
    ///格式化输出差分表与插值多项式
    const void Display() const;
    /// 获得存储的插值多项式
    const Polynomial<Real> getInterPoly()const{return interPoly_;};
    ///获得存储的差分表
    const TriangleTable<Real> getTableOfDividedDiffs()const{return tableOfDividedDiffs_;};
};

void NewtonInterp::QuickInterp()
{
    tableOfDividedDiffs_.clearAll();
    const InterpConditions IntpConds = this->getIntpConds();
    const RealBox IntpInfo = IntpConds.getIntpInfo();
    const RealVect RowNames = IntpConds.getIntpSets();
    tableOfDividedDiffs_.setRowNames(RowNames);
    const int RowNum = RowNames.size();
    tableOfDividedDiffs_.addEle(IntpInfo[0][1]);
    int SameSiteFlag = 1;
    int SiteIndex = 0;
    for (int tempRow = 1; tempRow < RowNum; tempRow++)
    {
        if (RowNames[tempRow - 1] == RowNames[tempRow])
        {
            SameSiteFlag++;
        }
        else
        {
            SameSiteFlag = 1;
            SiteIndex++;
        };
        for (int j = 0; j < SameSiteFlag; j++)
            tableOfDividedDiffs_.addEle(IntpInfo[SiteIndex][j + 1] / fac(j));
        for (int j = SameSiteFlag; j < tempRow + 1; j++)
            tableOfDividedDiffs_.addEle((tableOfDividedDiffs_.locate(tempRow, j - 1) - tableOfDividedDiffs_.locate(tempRow - 1, j - 1)) / (RowNames[tempRow] - RowNames[tempRow - j]));
    }
    this->GenerateNewtonPoly();
}

void NewtonInterp::AddIntpSiteThenInterp(const RealVect &_SiteInfo)
{
    const int oldRowNum = this->getIntpConds().getIntpSets().size();
    if (IntpConds_.SiteLegality(_SiteInfo))
    {
        IntpConds_.ResetIntpSite(_SiteInfo);
        const InterpConditions IntpConds = this->getIntpConds();
        const RealBox IntpInfo = IntpConds.getIntpInfo();
        const RealVect RowNames = IntpConds.getIntpSets();
        tableOfDividedDiffs_.setRowNames(RowNames);
        const int newRowNum = RowNames.size();
        int SameSiteFlag = 0;
        for (int tempRow = oldRowNum; tempRow < newRowNum; tempRow++)
        {
            SameSiteFlag++;
            for (int j = 0; j < SameSiteFlag; j++)
                tableOfDividedDiffs_.addEle(_SiteInfo[j + 1] / fac(j));
            for (int j = SameSiteFlag; j < tempRow + 1; j++)
                tableOfDividedDiffs_.addEle((tableOfDividedDiffs_.locate(tempRow, j - 1) - tableOfDividedDiffs_.locate(tempRow - 1, j - 1)) / (RowNames[tempRow] - RowNames[tempRow - j]));
        }
        this->MergeNewtonPoly(newRowNum - oldRowNum);
    }
    else
    {
        std::cout << "Warning: ILLEGAL INPUT." << std::endl;
    };
}

void NewtonInterp::GenerateNewtonPoly()
{
    const RealVect eleDiag = tableOfDividedDiffs_.getDiag();
    const RealVect RowNames = tableOfDividedDiffs_.getRowNames();
    const int RowNum = RowNames.size();
    Polynomial<Real> NewtonPolyTerm({1});
    Polynomial<Real> NewtonPoly({eleDiag[0]});
    for (int i = 1; i < RowNum; i++)
    {
        Polynomial<Real> tempPoly({-RowNames[i - 1], 1});
        NewtonPolyTerm = NewtonPolyTerm * tempPoly;
        NewtonPoly = NewtonPoly + eleDiag[i] * NewtonPolyTerm;
    }
    interPoly_ = NewtonPoly;
}

void NewtonInterp::MergeNewtonPoly(const int &_newAddedRowNum)
{
    const RealVect eleDiag = tableOfDividedDiffs_.getDiag();
    const RealVect RowNames = tableOfDividedDiffs_.getRowNames();
    const int RowNum = RowNames.size();
    Polynomial<Real> NewtonPolyTerm({1});
    Polynomial<Real> NewtonPoly({eleDiag[0]});
    for (int i = 1; i < RowNum - _newAddedRowNum; i++)
    {
        Polynomial<Real> tempPoly({-RowNames[i - 1], 1});
        NewtonPolyTerm = NewtonPolyTerm * tempPoly;
    }
    for (int i = RowNum - _newAddedRowNum; i < RowNum; i++)
    {
        Polynomial<Real> tempPoly({-RowNames[i - 1], 1});
        NewtonPolyTerm = NewtonPolyTerm * tempPoly;
        interPoly_ = interPoly_ + eleDiag[i] * NewtonPolyTerm;
    }
}

const void NewtonInterp::Display() const
{
    std::cout << tableOfDividedDiffs_;
    std::cout << interPoly_<<std::endl;
}

Real NewtonInterp::NevilleAitken(const Real &_x, const InterpConditions &_IntpConds)
{
    if (_IntpConds.is_PrimaryIntp())
    {
        const RealBox IntpInfo = _IntpConds.getIntpInfo();
        const int RowNum = IntpInfo.size();
        RealVect RowNames;
        RealVect tempCol;
        for (int i = 0; i < RowNum; i++)
        {
            RowNames.push_back(_x - IntpInfo[i][0]);
            tempCol.push_back(IntpInfo[i][1]);
        }
        for (int ColIndex = 1; ColIndex < RowNum; ColIndex++)
        {
            for (int RowIndex = 0; RowIndex < (RowNum - ColIndex); RowIndex++)
            {
                tempCol[RowIndex] = (RowNames[RowIndex] * tempCol[RowIndex + 1] - RowNames[RowIndex + ColIndex] * tempCol[RowIndex]) / (RowNames[RowIndex + ColIndex] - RowNames[RowIndex]);
            }
        }
        return tempCol[0];
    }
    else
    {
        std::cout << "Warning: ILLEGAL INPUT." << std::endl;
        return -1;
    }
};
#else
//do nothing
#endif